{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "success! :)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np \n",
    "import pandas as pd\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sksurv.datasets import get_x_y \n",
    "from sksurv.ensemble import GradientBoostingSurvivalAnalysis\n",
    "\n",
    "from sksurv.metrics import concordance_index_censored, integrated_brier_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#TODO 读取包含指标数据和融合数据的文件----------------------------------------------------------\n",
    "indicatorDataPath = './data/trainset.xlsx'\n",
    "indicatorData = pd.read_excel(indicatorDataPath)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "indicatorData['state'] = indicatorData['state'].astype(int).astype(bool)\n",
    "indicatorData['lifetime'] = indicatorData['lifetime'].astype(int) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#TODO 指标数据列名\n",
    "indicatorColumns = ['novelty','number_of_patents_cited', \n",
    "        'science_linkage','number_of_IPC4','number_of_IPC8',         \n",
    "       'number_of_claims', 'number_of_independent_claims','number_of_dependent_claims', \n",
    "       'number_of_priorities', 'review_duration',     \n",
    "       'number_of_patent_families_dwpi','number_of_patent_families_inpadoc',\n",
    "       'number_of_countries_applying',       \n",
    "       'number_of_inventors', 'number_of_patentees',\n",
    "       'number_of_inventor_patents', 'number_of_first_inventor_patents',\n",
    "       'number_of_patentee_patents', 'patentee_category',  \n",
    "       \n",
    "       'state','lifetime']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#TODO 数据集划分：训练集和验证集2000-2020年专利数据（预测集2021，此案例未展示）\n",
    "trainTestData = indicatorData[indicatorData['public_year']<2021]\n",
    "print(trainTestData.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "trainData, testData = train_test_split(trainTestData,test_size=0.2, random_state=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#TODO 构建生存预测模型数据集\n",
    "trainColumns = indicatorColumns.copy()  # 选用指标数据来训练生存预测模型\n",
    "trainTestSet = trainTestData.loc[:, trainColumns]\n",
    "trainSet = trainData.loc[:, trainColumns]\n",
    "testSet = testData.loc[:, trainColumns]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature, label = get_x_y(data_frame = trainTestSet, \n",
    "               attr_labels = ['state', 'lifetime'],\n",
    "               pos_label = True,  #event发生对应的值（专利失效-True）\n",
    "               survival=True)\n",
    "xTrain, yTrain = get_x_y(data_frame = trainSet, \n",
    "               attr_labels = ['state', 'lifetime'],\n",
    "               pos_label = True, \n",
    "               survival=True)\n",
    "xTest, yTest = get_x_y(data_frame = testSet, \n",
    "               attr_labels = ['state', 'lifetime'],\n",
    "               pos_label = True,  \n",
    "               survival=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>#sk-container-id-1 {color: black;background-color: white;}#sk-container-id-1 pre{padding: 0;}#sk-container-id-1 div.sk-toggleable {background-color: white;}#sk-container-id-1 label.sk-toggleable__label {cursor: pointer;display: block;width: 100%;margin-bottom: 0;padding: 0.3em;box-sizing: border-box;text-align: center;}#sk-container-id-1 label.sk-toggleable__label-arrow:before {content: \"▸\";float: left;margin-right: 0.25em;color: #696969;}#sk-container-id-1 label.sk-toggleable__label-arrow:hover:before {color: black;}#sk-container-id-1 div.sk-estimator:hover label.sk-toggleable__label-arrow:before {color: black;}#sk-container-id-1 div.sk-toggleable__content {max-height: 0;max-width: 0;overflow: hidden;text-align: left;background-color: #f0f8ff;}#sk-container-id-1 div.sk-toggleable__content pre {margin: 0.2em;color: black;border-radius: 0.25em;background-color: #f0f8ff;}#sk-container-id-1 input.sk-toggleable__control:checked~div.sk-toggleable__content {max-height: 200px;max-width: 100%;overflow: auto;}#sk-container-id-1 input.sk-toggleable__control:checked~label.sk-toggleable__label-arrow:before {content: \"▾\";}#sk-container-id-1 div.sk-estimator input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 div.sk-label input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 input.sk-hidden--visually {border: 0;clip: rect(1px 1px 1px 1px);clip: rect(1px, 1px, 1px, 1px);height: 1px;margin: -1px;overflow: hidden;padding: 0;position: absolute;width: 1px;}#sk-container-id-1 div.sk-estimator {font-family: monospace;background-color: #f0f8ff;border: 1px dotted black;border-radius: 0.25em;box-sizing: border-box;margin-bottom: 0.5em;}#sk-container-id-1 div.sk-estimator:hover {background-color: #d4ebff;}#sk-container-id-1 div.sk-parallel-item::after {content: \"\";width: 100%;border-bottom: 1px solid gray;flex-grow: 1;}#sk-container-id-1 div.sk-label:hover label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 div.sk-serial::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: 0;}#sk-container-id-1 div.sk-serial {display: flex;flex-direction: column;align-items: center;background-color: white;padding-right: 0.2em;padding-left: 0.2em;position: relative;}#sk-container-id-1 div.sk-item {position: relative;z-index: 1;}#sk-container-id-1 div.sk-parallel {display: flex;align-items: stretch;justify-content: center;background-color: white;position: relative;}#sk-container-id-1 div.sk-item::before, #sk-container-id-1 div.sk-parallel-item::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: -1;}#sk-container-id-1 div.sk-parallel-item {display: flex;flex-direction: column;z-index: 1;position: relative;background-color: white;}#sk-container-id-1 div.sk-parallel-item:first-child::after {align-self: flex-end;width: 50%;}#sk-container-id-1 div.sk-parallel-item:last-child::after {align-self: flex-start;width: 50%;}#sk-container-id-1 div.sk-parallel-item:only-child::after {width: 0;}#sk-container-id-1 div.sk-dashed-wrapped {border: 1px dashed gray;margin: 0 0.4em 0.5em 0.4em;box-sizing: border-box;padding-bottom: 0.4em;background-color: white;}#sk-container-id-1 div.sk-label label {font-family: monospace;font-weight: bold;display: inline-block;line-height: 1.2em;}#sk-container-id-1 div.sk-label-container {text-align: center;}#sk-container-id-1 div.sk-container {/* jupyter's `normalize.less` sets `[hidden] { display: none; }` but bootstrap.min.css set `[hidden] { display: none !important; }` so we also need the `!important` here to be able to override the default hidden behavior on the sphinx rendered scikit-learn.org. See: https://github.com/scikit-learn/scikit-learn/issues/21755 */display: inline-block !important;position: relative;}#sk-container-id-1 div.sk-text-repr-fallback {display: none;}</style><div id=\"sk-container-id-1\" class=\"sk-top-container\"><div class=\"sk-text-repr-fallback\"><pre>GradientBoostingSurvivalAnalysis(max_depth=5, n_estimators=140)</pre><b>In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. <br />On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.</b></div><div class=\"sk-container\" hidden><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-1\" type=\"checkbox\" checked><label for=\"sk-estimator-id-1\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">GradientBoostingSurvivalAnalysis</label><div class=\"sk-toggleable__content\"><pre>GradientBoostingSurvivalAnalysis(max_depth=5, n_estimators=140)</pre></div></div></div></div></div>"
      ],
      "text/plain": [
       "GradientBoostingSurvivalAnalysis(max_depth=5, n_estimators=140)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#TODO 构建与训练生存预测模型GBSA----------------------------------------------------------------\n",
    "surv_func = GradientBoostingSurvivalAnalysis(n_estimators=140,  # 树的数量\n",
    "                           loss=\"coxph\",\n",
    "                           learning_rate=0.1,\n",
    "                           max_depth=5)\n",
    "surv_func.fit(xTrain, yTrain)\n",
    "# 训练时间：17-30min"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "train concordance index: 0.7667344068695384\n",
      "test concordance index: 0.7333096980245946\n"
     ]
    }
   ],
   "source": [
    "# 模型评估\n",
    "# 一致性指数\n",
    "# train\n",
    "trainPre = surv_func.predict(xTrain)\n",
    "concordanceIndexTrain = concordance_index_censored(yTrain['state'],\n",
    "        yTrain['lifetime'], trainPre)\n",
    "print('train concordance index:' ,concordanceIndexTrain[0])\n",
    "# test\n",
    "TestPre = surv_func.predict(xTest)\n",
    "concordanceIndexTest = concordance_index_censored(yTest['state'],\n",
    "        yTest['lifetime'], TestPre)\n",
    "print('test concordance index:' ,concordanceIndexTest[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "range of time: lower: 28.0  ;upper: 235.0\n",
      "train integrated brier score: 0.1515796190231837\n",
      "test integrated brier score: 0.1607947584710229\n"
     ]
    }
   ],
   "source": [
    "# 综合Brier评分\n",
    "lower, upper = np.percentile(label['lifetime'], [10, 95])  #计算分位数\n",
    "print('range of time: lower:',lower,' ;upper:', upper)\n",
    "timesIBS = np.arange(lower, upper + 1)\n",
    "# train\n",
    "predsTrain = np.row_stack([fn(timesIBS) for fn in surv_func.predict_survival_function(xTrain)])\n",
    "trainIBS = integrated_brier_score(label, yTrain, predsTrain, timesIBS)\n",
    "print('train integrated brier score:', trainIBS)\n",
    "# test\n",
    "predsTest = np.row_stack([fn(timesIBS) for fn in surv_func.predict_survival_function(xTest)])\n",
    "testIBS = integrated_brier_score(label, yTest, predsTest, timesIBS)\n",
    "print('test integrated brier score:', testIBS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('index of 10 year in event time:', list(surv_func.event_times_).index(10*12))\n",
    "surv_func.event_times_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 2011\n"
     ]
    }
   ],
   "source": [
    "# TODO 确定生存风险函数-------------------------------------------------------------------\n",
    "# 高价值专利寿命阈值，及相应的公开年\n",
    "lifetimeSplit = 10  #维持时间超过10年判定为高价值专利\n",
    "publicYearSplit = 2021-lifetimeSplit  #此公开年之前的专利可用来进行阈值验证\n",
    "print(lifetimeSplit, publicYearSplit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of train data public before 2012: 4129\n",
      "lifetime over 10 year number: 3084\n"
     ]
    }
   ],
   "source": [
    "# 训练集中维持时间超过10年的占比\n",
    "trainDataCopy = trainData.copy()\n",
    "trainDataBefore2011 = trainDataCopy[trainDataCopy['public_year']<=publicYearSplit]\n",
    "print('number of train data public before 2012:',trainDataBefore2011.shape[0])\n",
    "lifetimeOver10year_number = trainDataBefore2011[trainDataBefore2011['lifetime']>lifetimeSplit*12].shape[0]\n",
    "print('lifetime over {} year number:'.format(lifetimeSplit), lifetimeOver10year_number) \n",
    "# 选定生存评先阈值为 测试集2000-2011年专利风险评分为升序排序74.69%对应的风险评分\n",
    "# 后面按照这个比例试一下就可以，也可以写代码搜索。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "#TODO 测试集（所有）预测结果（风险评分，第十年生存概率）和验证集原始数据保存到excel\n",
    "# 风险评分\n",
    "testRiskScore = surv_func.predict(xTest)\n",
    "# 第10年生存概率\n",
    "position10year = list(surv_func.event_times_).index(lifetimeSplit*12)\n",
    "survival_function_test = surv_func.predict_survival_function(xTest, return_array=True)\n",
    "testSurvivalProbability10year = survival_function_test[:,position10year]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "high value number: 4325\n",
      "other number: 1802\n"
     ]
    }
   ],
   "source": [
    "#TODO 设定风险评分阈值，划分是否为高价值专利\n",
    "higValueThreshold = 0.065\n",
    "testPreHighValueLabel = [1 if score<higValueThreshold else 0 \n",
    "                         for score in testRiskScore]\n",
    "print('high value number:', testPreHighValueLabel.count(1))\n",
    "print('other number:', testPreHighValueLabel.count(0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "testDataCopy = testData.copy()\n",
    "testDataCopy['risk_score'] = testRiskScore\n",
    "testDataCopy['survival_probability'] = testSurvivalProbability10year\n",
    "testDataCopy['high_value_pred'] = testPreHighValueLabel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of test data public before 2012: 1032\n",
      "lifetime over 10 year number: 767\n",
      "predict high value number: 771\n",
      "predict correct number: 597\n"
     ]
    }
   ],
   "source": [
    "#TODO 利用 测试集 中2000-2011年专利看划分的高价值专利阈值是否合理。\n",
    "# （预测为高价值与真实超过10年的两个集合，交集占比）\n",
    "testDataBefore2011 = testDataCopy[testDataCopy['public_year']<=publicYearSplit]\n",
    "print('number of test data public before {}:'.format(publicYearSplit+1),\n",
    "      testDataBefore2011.shape[0])\n",
    "# 维持时间超过10年专利数量\n",
    "lifetimeOver10year_number = testDataBefore2011[testDataBefore2011['lifetime']>lifetimeSplit*12].shape[0]\n",
    "# 被预测为高价值专利数量\n",
    "predHighValue_number = testDataBefore2011[testDataBefore2011['high_value_pred']==1].shape[0]\n",
    "# 被预测为高价值专利数量且维持时间超过10年\n",
    "predCorrect_number = testDataBefore2011[(testDataBefore2011['high_value_pred']==1)&(\n",
    "    testDataBefore2011['lifetime']>lifetimeSplit*12)].shape[0]\n",
    "\n",
    "print('lifetime over {} year number:'.format(lifetimeSplit), lifetimeOver10year_number)  # 验证集中符合条件的\n",
    "print('predict high value number:', predHighValue_number)\n",
    "print('predict correct number:', predCorrect_number)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py38",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "20a1dc4dfb8b30f987c459d716b3c8ae463117e9025402208219e2a6b3d7e9d2"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
